Casimir forces in the time domain: I. Theory 
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We introduce a method to compute Casimir forces in arbitrary geometries and for arbitrary 
materials based on the finite-difference time-domain (FDTD) scheme. The method involves the 
time-evolution of electric and magnetic fields in response to a set of current sources, in a modified 
medium with frequency-independent conductivity. The advantage of this approach is that it allows 
one to exploit existing FDTD software, without modification, to compute Casimir forces. In this 
manuscript, part I, we focus on the derivation, implementation choices, and essential properties of 
the time-domain algorithm, both considered analytically and illustrated in the simplest parallel-plate 
geometry. Part II presents results for more complex two- and three-dimensional geometries. 



I. INTRODUCTION 

In recent years, Casimir forces arising from quantum 
vacuum fluctuations of the electromagnetic field [TJ EJ |3] 
have become the focus of intense theoretical and ex- 
perimental effort nnEnzuHnnniniiiijiniinjiTi 

us Ha [m [H Ha Uni ini- This effect has been veri- 
fied via many experiments [22j [23j [Ml El], most com- 
monly in simple, one-dimensional geometries involving 
parallel plates or approximations thereof, with some ex- 
ceptions [26] . A particular topic of interest is the ge- 
ometry and material dependence of the force, a sub- 
ject that has only recently begun to be addressed in ex- 
periments |26j and by promising new theoretical meth- 
ods EZll2ll|23|3aini[3a[MllMl33IMll3Zl|38]. 
For example, recent works have shown that it is possi- 
ble to find unusual effects arising from many-body in- 
teractions or from systems exhibiting strongly coupled 
material and geometric dispersion [39j HQI [HI |42l |43] . 
These numerical studies have been mainly focused in two- 
dimensional ^13. ^44, 45, 46 or simple three-dimensional 
constant-cross-section geometries [331 HOI HZ] for which 
numerical calculations are tractable. 

In this manuscript, we present a simple and general 
method to compute Casimir forces in arbitrary geome- 
tries and for arbitrary materials that is based on a 
finite-difference time-domain (FDTD) scheme in which 
Maxwell's equations are evolved in time [48]. A time- 
domain approach offers a number of advantages over 
previous methods. First, and foremost, it enables re- 
searchers to exploit powerful free and commercial FDTD 
software with no modification. The generality of many 
available FDTD solvers provides yet another means to 
explore the material and geometry dependence of the 
force, including calculations involving anisotropic di- 
electrics [49] and/or three-dimensional problems. Sec- 
ond, this formulation also offers a fundamentally different 
viewpoint on Casimir phenomena, and thus new oppor- 
tunities for the theoretical and numerical understanding 
of the force in complex geometries. 

Our time-domain method is based on a standard for- 
mulation in which the Casimir force is expressed as 
a contour integral of the frequency-domain stress ten- 



sor [2]. Like most other methods for Casimir calculations, 
the stress tensor method typically involves evaluation at 
imaginary frequencies, which we show to be unsuitable 
for FDTD. We overcome this difficulty by exploiting a 
recently-developed exact equivalence between the system 
for which we wish to compute the Casimir force and a 
transformed problem in which all material properties are 
modified to include dissipation I^T)!. To illustrate this 
approach, we consider a simple choice of contour, corre- 
sponding to a conductive medium, that leads to a simple 
and efficient time-domain implementation. Finally, using 
a free, widely-available FDTD code [5T], we compute the 
force between two vacuum-separated perfectly-metallic 
plates, a geometry that is amenable to analytical calcu- 
lations and which we use to analyze various important 
features of our method. An illustration of the power and 
fiexibility of this method will be provided in a subse- 
quent article |52j . currently in preparation, in which we 
will demonstrate computations of the force in a number 
of non-trivial (dispersive, three-dimensional) geometries 
as well as further refinements to the method. 



II. METHOD 

In what follows, we derive a numerical method to 
compute the Casimir force on a body using the FDTD 
method. The basic steps involved in computing the force 
are: 

(1) Map the problem exactly onto a new problem with 
dissipation given by a frequency-independent con- 
ductivity a. 

(2) Measure the electric E and magnetic H fields in 
response to current pulses placed separately at each 
point along a surface enclosing the body of interest. 

(3) Integrate these fields in space over the enclosing 
surface and then integrate this result, multiplied by 
a known function g(^t), over time t, via Eq. (29 1. 

The result of this process is the exact Casimir force 
(in the limit of sufficient computational resolution), ex- 
pressed via Eq. ( 29 I and requiring only the time-evolution 

of Eqs. (ITsHTel). 
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In this section, we describe the mathematical develop- 
ment of our time-domain computational method, start- 
ing from a standard formulation in which the Casimir 
force is expressed as a contour integral of the frequency- 
domain stress tensor. We consider the frequency do- 
main for derivation purposes only, since the final tech- 
nique outlined above resides entirely in the time domain. 
In this framework, computing the Casimir force involves 
the repeated evaluation of the photon Green's function 
Gij over a surface S surrounding the object of interest. 
Our goal is then to compute Gij via the FDTD method. 
The straightforward way to achieve this involves com- 
puting the Fourier transform of the electric field in re- 
sponse to a short pulse. However, in most methods a 
crucial step for evaluating the resulting frequency inte- 
gral is the passage to imaginary frequencies, correspond- 
ing to imaginary time. We show that, in the FDTD 
this, gives rise to exponentially growing solutions and is 
therefore unsuitable. Instead, we describe an alternative 
formulation of the problem that exploits a recently pro- 
posed equivalence in which contour deformations in the 
complex frequency-domain w(^) correspond to introduc- 
ing an effective dispersive, dissipative medium at a real 
"frequency" ^. From this perspective, it becomes simple 
to modify the FDTD Maxwell's equations for the pur- 
pose of obtaining well-behaved stress tensor frequency 
integrands. We illustrate our approach by considering 
a contour corresponding to a medium with frequency- 
independent conductivity a. This contour has the ad- 
vantage of being easily implemented in the FDTD, and 
in fact is already incorporated in most FDTD solvers. 
Finally, we show that it is possible to abandon the fre- 
quency domain entirely in favor of evaluating the force 
integral directly in the time domain, which offers several 
conceptual and numerical advantages. 



and {H,{y,lu)Hj{t\uj)): 



+ e(r, f {E,{r) E,{v))^ ~ ^ {E^,{v) E,{r)) 



(2) 



where both the electric and magnetic field correlation 
functions can be written as derivatives of a vector poten- 
tial operator A'^(r,a;): 

E,{r,uj) = -iujAf{r,uj) 



(3) 
(4) 



We explicitly place a superscript on the vector poten- 
tial in order to refer to our choice of gauge [Eqs. ([s]- 
|4])], in which E is obtained as a time-derivative of A. 
The fluctuation-dissipation theorem relates the corre- 
lation function of A^ to the photon Green's function 



{Af{r,u;)Af{r',Lu)) = Im G^(c., r, r'), (5) 

where Gfj is the vector potential Af in response to an 
electric dipole current J along the direction: 

1 



V X — X -uj'^e(r,uj) 

/i(r,w) 



Gf(c^;r,r')-(5(r~r')e,, 
(6) 

Given Gfj, one can use Eqs. P|j4| in conjunction with 
Eq. (|5| to express the field correlafion functions at points 
r and r' in terms of the photon Green's function: 



{E,{r,uj)E,{r',u;)) - -u;^ lmGfjLj;r,r') 



(7) 



{H,{r,u;)H,{r',u;)) = — (Vx),KV'x),„,ImGf„(r,r',c^), 

TT 

(8) 



Stress Tensor Formulation 



The Casimir force on a body can be expressed [5] as an 
integral over any closed surface S (enclosing the body) of 
the mean electromagnetic stress tensor {Tij{r,uj)). Here 
r denotes spatial position and u; frequency. In particular, 
the force in the ith direction is given by: 

poo 

F,^ dcuffY.('^,,{r,cu))dS,, (1) 



The stress tensor is expressed in terms of correlation 
functions of the the field operators {Ei{r,uj)Ej{r',uj)) 



In order to find the force via Eq. ([T]), we must first 
compute Gfj{r,r' = r,w) at every r on the surface of 
integration S, and for every uj [5]. Equation ^ can 
be solved numerically in a number of ways, such as by 
a finite-difference discretization [30]: this involves dis- 
cretizing space and solving the resulting matrix eigen- 
value equation using standard numerical linear algebra 
techniques [53] |54| • We note that finite spatial discretiza- 
tion automatically regularizes the singularity in Gf^ at 
r = r', making Gfj finite everywhere [30j . 
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B. Complex Frequency Domain 



C. Time Domain Approach 



The present form of Eq. ^ is of limited computational 
utility because it gives rise to an oscillatory integrand 
with non-negligible contributions at all frequencies, mak- 
ing numerical integration difficult |30j . However, the in- 
tegral over Lo can be re-expressed as the imaginary part 
of a contour integral of an analytic function by commut- 
ing the uj integration with the Im operator in Eqs. ([^jsj. 
Physical causality implies that there can be no poles in 
the integrand in the upper complex plane. The inte- 
gral, considered as a complex contour integral, is then 
invariant if the contour of integration is deformed above 
the real frequency axis and into the first quadrant of the 
complex frequency plane, via some mapping uj — > ^(0- 
This allows us to add a positive imaginary component 
to the frequency, which causes the force integrand to de- 
cay rapidly with increasing ^ |50j . In particular, upon 
deformation, Eq. ([6| is mapped to: 



V X 



1 



-Vx -a>2(^)e(r,u;) 



Gf(e;r,r')-<5(r-r')e,: 
(9) 



/i(r,w) 

and Eqs. (|7||8]) are mapped to: 



{H,{r,u;)H,{r',u;)) = — (Vx),a(V'x),„GL(r,r',u;), 

TT 



(10) 
(11) 



Equation ([T]) becomes: 

i^.^Imy d^^ § Y.(T,,{r,co))dS,, (12) 

^ surface j 

[Note that a finite spatial grid (as used in the present 
approach) requires no further regularization of the inte- 
grand, and the finite value of all quantities means there 
is no difficulty in commuting the Im operator with the 
integration.] 

We can choose from a general class of contours, pro- 
vided that they satisfy w(0) = and remain above 
the real ^ axis. The standard contour w(^) = is a 
Wick rotation, which is known to yield a force integrand 
that is smooth and exponentially decaying in ^ [5]. In 
general, the most suitable contour will depend on the 
numerical method being employed. A Wick rotation 
guarantees a strictly positive-definite and real-symmetric 
Green's function, making Eq. ^ solvable by the most ef- 
ficient numerical techniques (e.g. the conjugate-gradient 
method) [54j- One can also solve Eq. ([6| for arbitrary 
^(0 [Sn], but this will generally involve the use of di- 
rect solvers or more complicated iterative techniques [55] . 
However, the class of contours amenable to an efficient 
time-domain solution is more restricted. For instance, a 
Wick rotation turns out to be unstable in the time do- 
main because it implies the presence of gain |50j . 



It is possible to solve Eq. (|6| in the time domain 
by evolving Maxwell's equations in response to a delta- 
function current impulse J(r,t) = 5{r — r')S{t — t')ej in 
the direction of e^ . Gfj can then be directly computed 
from the Fourier transform of the resulting E field. How- 
ever, obtaining a smooth and decaying force integrand 
requires expressing the mapping lu — > a;(^) in the time- 
domain equations of motion. A simple way to see the 
effect of this mapping is to notice that Eq. ^ can be 
viewed as the Green's function at real "frequency" ^ and 
complex dielectric |5U] : 



\0 



e(r) 



(13) 



where for simplicity we have taken /i and e to be 
frequency-independent. We assume this to be the case 
for the remainder of the manuscript. At this point, it is 
important to emphasize that the original physical system 
£ at a frequency lo is the one in which Casimir forces and 
fluctuations appear; the dissipative system Sc at a fre- 
quency ^ is merely an artificial technique introduced to 
compute the Green's function. 

Integrating along a frequency contour w(^) is therefore 
equivalent to making the medium dispersive in the form 
of Eq. (13 1. Consequently, the time domain equations 



of motion under this mapping correspond to evolution 
of the fields in an effective dispersive medium given by 
ec(r,C). 

To be suitable for FDTD, this medium should have 
three properties: it must respect causality, it cannot sup- 
port gain (which leads to exponential blowup in time- 
domain), and it should be easy to implement. A Wick 
rotation is very easy to implement in the time-domain: it 
corresponds to setting Sc = — £. However, a negative ep- 
silon represents gain (the refractive index is i-y/e, where 
one of the signs corresponds to an exponentially grow- 
ing solution) . We are therefore forced to consider a more 
general, frequency-dependent Sc- 

Implementing arbitrary dispersion in FDTD generally 
requires the introduction of auxiliary fields or higher or- 
der time-derivative terms into Maxwell's equations, and 
can in general become computationally expensive |48| . 
The precise implementation will depend strongly on the 
choice of contour a;(^). However, almost any dispersion 
will suit our needs, as long as it is causal and dissipative 
(excluding gain). A simple choice is an £c(r,^) corre- 
sponding to a medium with frequency-independent con- 
ductivity cr: 



£e(r,0=e(r) 1 



J 



(14) 



This has three main advantages: first, it is imple- 
mented in many FDTD solvers currently in use; second, 
it is numerically stable; and third, it can be efficiently im- 
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plemented without an auxiliary differential equation |48j . 
In this case, the equations of motion in the time domain 
are given by: 



cussed in the Appendix, this choice of vector potential 
implies a frequency-independent magnetic conductivity, 
and a magnetic, instead of electric, current. The result- 
ing time-domain equations of motion are: 



dt 



-V X E 

V X H - creE - J 



(15) 
(16) 



9eE 



Writing the conductivity term as ae is slightly nonstan- 
dard, but is convenient here for numerical reasons. In 
conjunction with Eqs. ([3jj4]), and a Fourier transform in 
^, this yields a photon Green's function given by: 



dt 



-V X E + (T/iH - J (20) 
V X H (21) 



In this gauge, the new photon Green's function = 

(Af (r,0^f (r',0) and the field H in response to the 
current source J are related by: 



Vx^tVx ~ee{r) ( 1 + ^ 



G,(e;r,r') = 5(r-r')e, 
(17) 



This corresponds to picking a frequency contour of the 
form: 



(18) 



Note that, in the time domain, the frequency of the fields 
is ^, and not to, i.e. their time dependence is e~'^*. The 
only role of the conductivity a here is to introduce an 



imaginary component to Eq. ( 17 1 in correspondence with 



a complex-frequency mapping. It also explicitly appears 
in the final expression for the force, Eq. ( 12 1, as a multi- 
plicative (Jacobian) factor. 

The standard FDTD method involves a discretized 



form of Eqs. (15-161, from which one obtains E and 
B, not Gfy However, in the frequency domain, the 
photon Green's function, being the solution to Eq. 
solves exactly the same equations as those satisfied by the 
electric field E, except for a simple multiplicative factor 
in Eq. (jsj). Specifically, Gfj is given in terms of E by: 



(19) 



where Eij{r,^) denotes the field in the ith direction due 
to a dipole current source J(r, t) ~ J'{t)6{r — r')cj placed 
at r' with time-dependence J^{t), e.g. J^{t) = d{t). 

In principle, we can now compute the electric- and 



magnetic-field correlation functions by using Eqs. (10 
11 1, with cj(^) given by Eq. ( 18 1, and by setting r — rin 
Eq. (111. Since we assume a discrete spatial grid, no sin- 
gularities arise for r = r', and in fact any r-independent 
contribution is canceled upon integration over S. This is 
straightforward for Eq. ([t]), since the E-field correlation 
function only involves a simple multiplication by uj'^{£,)- 
However, the H-field correlation function, Eq. ([s]), in- 
volves derivatives in space. Although it is possible to 
compute these derivatives numerically as finite differ- 
ences, it is conceptually much simpler to pick a differ- 
ent vector potential, analogous to Eqs. ([3}]4|, in which H 
is the time-derivative of a vector potential A-^. As dis- 



where the magnetic-field correlation function: 

{H,ir,OH,ir',0) = -^2(0G^(e r, r'), (23) 

TT 

is now defined as a frequency multiple of rather than 
by a spatial derivative of Gfj . 

This approach to computing the magnetic correlation 
function has the advantage of treating the electric and 
magnetic fields on the same footing, and also allows us 
to examine only the field response at the location of the 
current source. The removal of spatial derivatives also 
greatly simplifies the incorporation of discretization into 
our equations (see Appendix for further discussion) . The 
use of magnetic currents and conductivities, while un- 
physical, are easily implemented numerically. Alterna- 
tively, one could simply interchange e and E and H, 
and run the simulation entirely as in Eqs. ( 15 16 1. 



The full force integral is then expressed in the sym- 
metric form: 



rfC5(0(rf(0 + rf(0), (24) 



where 



rf(0 ^ #E=W (i?.,(r)-^%E^^.'=w) « 

S j \ k ) 



rf(0 



represent the surface-integrated field responses in the fre- 
quency domain, with Eij{r) = Eij{r;£^). For notational 
simplicity, we have also defined: 



m 



0(0 



(27) 



Here, the path of integration has been extended to the 
entire real ^-axis with the use of the unit-step function 
0(^) for later convenience. 
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The product of the fields with g{^) naturally decom- 
poses the problem into two parts: computation of the 
surface integral of the field correlations T, and of the 
function .g(^). The Ti contain all the structural informa- 
tion, and are straightforward to compute as the output 
of any available FDTD solver with no modification to the 
code. This output is then combined with g{£,), which is 
easily computed analytically, and integrated in Eq. (24) 
to obtain the Casimir force. As discussed in Sec. IV A 
the effect of spatial and temporal discretization enters ex- 
plicitly only as a slight modification to g{£^) in Eq. (24 1 
leaving the basic conclusions unchanged. 



D. Evaluation in the Time Domain 



It is straightforward to evaluate Eq. (24 1 in the fre 



quency domain via a dipole current J^{t) — S(t), which 
yields a constant-amplitude current J^{S,) = 1. Using the 
frequency-independent conductivity contour Eq. (18 1, 



corresponding to Eqs. (15 
plicit form for (?(^): 



16 1, we find the following ex- 



J 



1 + 



0(0 (28) 



One important feature of Eq. (28 1 is that g(^) ^Ao^VC 
becomes singular in the limit as ^ ^ 0. Assuming that 
r^(^) and r^(^) are continuous at ^ = (in general they 
will not be zero), this singularity is integrable. However, 
it is cumbersome to integrate in the frequency domain, 
as it requires careful consideration of the time window 
for calculation of the field Fourier transforms to ensure 
accurate integration over the singularity. 

As a simple alternative, we use the convolution theo- 
rem to re-express the frequency (^) in tegr al of the prod- 
uct oi g{£,) and r-^(^) arising in Eq. (24) as an integral 
over time t of their Fourier transform&^~g(— i) and T^{t). 
Technically, the Fourier transform of g{£,) does not exist 
because g{(,) ~ f for large ^. However, the integral is 
regularized below using the time discretization, just as 
the Green's function above was regularized by the spa- 
tial discretization. (As a convenient abuse of notation, ^ 
arguments will always denote functions in the frequency 
domain, and t arguments their Fourier transforms in the 
time domain.) 

Taking advantage of the causality conditions 
{T^{t), T"{t) = for t < 0) yields the following 
expression for the force expressed purely in the time 
domain: 



F, = Im- 



dt gi-t) (Tfit) + Tf{t)) 



(29) 



The advantage of evaluating the force integral in the 
time domain is that, due to the finite conductivity and 
lack of sources for t > 0, T{t) will rapidly decay in time. 
As will be shown in the next section, g{—t) also decays 



with time. Hence, although dissipation was originally 
introduced to permit a natural high-frequency cutoff to 
our computations, it also allows for a natural time cutoff 
T. We pick T such that, for times t > T, knowledge 
of the fields will not change the force result in Eq. (29) 



beyond a predetermined error threshold. This approach 
is very general as it requires no precise knowledge of how 
the fields decay with time. 



E. Properties of g{^t) 

Given g{S,), the desired function g{—t) is a Fourier 
transform. However, the discretization of time in FDTD 
implies that the frequency domain becomes periodic and 
that g{t) = g{n/S.t) are actually Fourier series coefficients, 
given by: 



g{nAt)^ / d^gdiOe 
Jo 



(30) 



where gd{0 is the discretized form of Eq. (27) and is 



given in the Appendix by Eq. ( 38 1 . These Fourier series 



coefficients are computed by a sequence of numeric in- 
tegrals that can be evaluated in a variety of ways. It is 
important to evaluate them accurately in order to resolve 
the effect of the ^ = singularity. For example, one could 
use a Clenshaw-Curtis scheme developed specifically for 
Fourier integrals [55' , or simply a trapezoidal rule with a 
large number of points that can be evaluated relatively 
quickly by an FFT (e.g. for this particular g{^), 10'' 
points is sufficient). 

Since it is possible to employ strictly-real current 
sources in FDTD, giving rise to real F, and since we 
are only interested in analyzing the influence of g(t) on 
Eq. (29 1, it suffices to look at Im g{~t). Furthermore, g{t) 
will exhibit rapid oscillations at the Nyquist frequency 
due to the delta- function current, and therefore it is more 
convenient to look at its absolute value. Figure [T] below, 
plots the envelope of | Im g{—t)\ as a function of t, where 
again, g{t) is the Fourier transform of Eq. ( p7| . 

As anticipated in the previous section, g{t) decays in 
time. Interestingly, it exhibits a transition from ^ 
decay at ti = to ~ i^^/^ decay for large a. The slower 
decay at long times for larger a arises from a transition 
in the behavior of Eq. (28 1 from the singularity at ^ = 0. 



III. PROPERTIES OF THE METHOD 

In this section we discuss the practical implementation 
of the time-domain algorithm (using a freely-available 
time domain solver [51 that required no modification). 
We analyze its properties applied to the simplest parallel- 
plate geometry [Fig. [2] , which illustrate the essential fea- 
tures in the simplest possible context. In particular, we 
analyze important computational properties such as the 
convergence rate and the impact of different conductiv- 




FIG. 1: I lmg{t)\ for various values of a, illustrating the tran- 
sition from to t~^^^ power-law decay as a increases. Be- 
cause there are strong oscillations in g{t) at the Nyquist fre- 
quency for intermediate a, for clarity we plot the positive and 
negative terms in g(t) as separate components. 



ity choices. Part II of this manuscript, in preparation, 
demonstrates the method for more comphcated two- and 
three-dimensional geometries [52 . 



A. Fields in Real Time 

The dissipation due to positive cr implies that the 
fields, and hence T^{t), will decay exponentially with 
time. Below, we use a simple one-dimensional exam- 
ple to understand the consequences of this dissipation 
for both the one-dimensional parallel plates and the two- 
dimensional piston configuration. The simplicity of the 
parallel-plate configuration allows us to examine much 
of the behavior of the time-domain response analytically. 
(The understanding gained from the one-dimensional ge- 
ometry can be applied to higher dimensions.) Further- 
more, we confirm that the error in the Casimir force due 
to truncating the simulation at finite time decreases ex- 
ponentially (rather than as t~^, as it would for no dissi- 
pation) . 



1. One- dimensional Parallel Plates 

To gain a general understanding of the behavior of the 
system in the time domain, we first examine a simple 
configuration of perfectly metallic parallel plates in one 
dimension. The plates are separated by a distance h (in 
units of an arbitrary distance a) in the x dimension, as 
shown by the inset of Fig. [2] The figure plots the field 
response T^{t) + T^{t), in arbitrary units, to a current 
source J^{t) = d{t) for increasing values of h, with the 
conductivity set at cr = 10 (27rc/a) . 

Figure[2]shows the general trend of the field response as 
a function of separation. For short times, all fields follow 




\ h=3 h=5 



20 40 60 80 100 120 140 160 180 200 

t(a/c) 

FIG. 2: r^(t) + r^(t) for a set of one-dimensional parallel 
plates as the separation h is varied. The inset shows the 
physical setup. 



the same power-law envelope, and later rapidly transition 
to exponential decay. Also plotted for reference is a 
curve, demonstrating that the envelope is in fact a power 
law. 

We can understand the power law envelope by con- 
sidering the vacuum Green's function in the case 
h ^ oo (analogous conclusions hold for G^). In the 
case h — > oo, one can easily solve for the vacuum Green's 
function G^(^, r — r') in one dimension for real frequency 



G^(e,r-r') 



ii\r-r'\ 



(31) 



We then analytically continue this expression to the 



complex frequency domain via Eq. ( 18 1 and compute the 
Fourier transform J d^e^^^G^ {oj {£,)). Setting r = r' in the 
final expression, one finds that, to leading order, G^{t) ^ 
^-3/2^ This explains the behavior of the envelope in Fig. [2] 
and the short-time behavior of the Green's functions: it 
is the field response of vacuum. 

Intuitively, the envelope decays only as a power in 
t because it receives contributions from a continuum 
of modes, all of which are individually decaying expo- 
nentially (this is similar to the case of the decay of 
correlations in a thermodynamic system near a critical 
point [.56J). For a finite cavity, the mode spectrum is 
discrete — the poles in the Green's function of the non- 
dissipative physical system are pushed below the real fre- 
quency axis in this dissipative, unphysical system, but 
they remain discretely spaced. 

At short times, the field response of a finite cav- 
ity will mirror that of an infinite cavity because the 
fields have not yet propagated to the cavity walls and 
back. As t increases, the cavity response will transi- 
tion to a discrete sum of exponentially decaying modes. 
From Eq. (18 1, higher- frequency modes have a greater 



imaginary-frequency component, so at sufficiently long 
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FIG. 3: Partial force as defined in Eq. (32 1 for one 



dimensional parallel plates as a function of time t. (Inset): 
Relative error A{t) as a function of t on a semi log scale. 



times the response will decay exponentially, the decay 
being determined by the lowest-frequency cavity mode. 
The higher the frequency of that mode, the faster the 
dissipation. 

This prediction is confirmed in Fig. |2] as /i decreases, 
the source "sees" the walls sooner. From the standpoint 
of computational efficiency, this method then works best 
when objects are in close proximity to one another (al- 
though not so close that spatial resolution becomes an 
issue), a situation of experimental interest. 



2. Convergence of the Force 

We now examine the force on the parallel plates. From 
the above discussions of the field decay and the decay of 
g{t), we expect the time integral in Eq. (29 1 to eventually 
converge exponentially as a function of time. In the in- 
terest of quantifying this convergence, we define the time 
dependent "partial force" Fi{t) as: 



F,;(t) = Im- 



dt' g{-t'){Tf{t') + Tfit')) (32) 



Letting Fi{oo) denote the t — > oo limit of Fi{t), which 
is the actual Casimir force, we define the relative error 
Ai{t) in the i-th component of the force as: 



F,it)-Fi{<x) 



Fiioo) 



(33) 



We plot F^{t) in Fig.[3]for the one-dimensional parallel- 
plate structure with different values of a. The inset plots 
A{t) for the same configuration. As expected, the asymp- 
totic value of Fx{t) is independent of a, and A(<) con- 
verges exponentially to zero. 

For a near zero, the force is highly oscillatory. In 
one dimension this gives the most rapid convergence 



with time, but it is problematic in higher dimensions. 
This is because, in higher-dimensional systems, S con- 
sists of many points, each contributing a response term 
as in Fig.|3] If a is small, every one of these terms will be 
highly oscillatory, and the correct force Eq. ( 32 1 will only 



be obtained through delicate cancellations at all points 
on S. Small cr is thus very sensitive to numerical error. 

Increasing a smooths out the response functions, as 
higher frequency modes are damped out. However, some- 
what counterintuitively, it also has the effect of slowing 
down the exponential convergence. One can understand 
the asymptotic behavior of the force by considering the 
equations of motion Eq. ( 17 1 as a function of a and ^. 
When the response function exhibits few if any oscilla- 
tions we are in the regime where a 3> ^. In this limit, 
the approximate equations of motion are: 



V X — 7-T-V X — iaS,e{r) 



M(r) 

In the limit of Eq. pil) 



G,(C;r,r') = <5(r~r')e, 



(34) 

the eigenfrequency ^ of a 
given spatial mode scales proportional to —ija. The 
lowest-frequency mode therefore has a time-dependence 
~ e"'-^*/'^, for some constant C > 0. Since the decay of 
the force at long times is determined by this mode, we ex- 
pect the decay time to scale inversely with u in the limit 
of very high tr. This is suggested in Fig. |3] and confirmed 
by further numerical experiments. 



Additionally, from Eq. ( 34 1 we see that in the case of a 



homogeneous one-dimensional cavity, the solutions have 
a quadratic dispersion ^ ~ il? , for spatial dependence 
e*'^^, and so the lowest cavity frequency scales as the in- 
verse square of the cavity size. This means that the rate 
of exponential convergence ofFig.[2]should vary as ~ ft 
in the limit of very large a. This scaling is approximately 
apparent from Fig. [2] and further experiments for much 
larger tr confirm the scaling. We thus see that in this 
limit, the effect of increasing a by some factor is analo- 
gous to increasing the wall spacing of the cavity by the 
square root of that factor. 

The present analysis shows that there are two unde- 
sirable extremes. When a is small, rapid oscillations in 
Fi(i)) will lead to large numerical errors in more than one 
dimension. When a is large, the resulting frequency shift 
will cause the cavity mode to decay more slowly, result- 
ing in a longer run time. The optimal a lies somewhere 
in between these two extremes and will generally depend 
on the system being studied. For the systems considered 
in this paper, with a typical scale « a, cr ^ 1 (27rc/a) 
appears to be a good value for efficient and stable time- 
domain computation. 



IV. CONCLUDING REMARKS 

An algorithm to compute Casimir forces in FDTD 
has several practical advantages. FDTD algorithms that 
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solve Maxwell's equations with frequency-independent 
conductivity, and even more complicated dispersions, are 
plentiful and well-studied. They are stable, convergent, 
and easily parallelized. Although the current formulation 
of our method requires the evaluation of Gij{v) along a 
surface S, requiring a separate calculation of the fields for 
each dipole source in S, all of these sources can be sim- 
ulated in parallel, with no communication between dif- 
ferent simulations until the very end of the computation. 
In addition, many FDTD solvers will allow the compu- 
tational cell for each source to be parallelized, providing 
a powerful method capable of performing large computa- 
tions. 
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APPENDIX 
A. Effects of Discretization 



The calculations of this paper employed non-dispersive 
materials in the original (w) system. However, the theo- 
retical analysis applies equally well to materials of arbi- 
trary dispersion. Any materials that can be implemented 
in an FDTD solver (e.g. a sum of Lorentzian dielec- 
tric resonances ^JS]) can also be included, and existing 
algorithms have demonstrated the ability to model real 
materials |1H1[S7]. Existing FDTD implementations also 
handle anisotropy in e and /i, multiple types of boundary 
conditions, and other complications [H]. 

In principle, the computational scaling of this FDTD 
method is comparable to finite-difference frequency- 
domain (FDFD) methods [30*. In both cases, each solver 
step (either a time step for FDTD or an iterative-solver 
step for FDFD) requires 0{N) work for N grid points. 
The number of time steps required by an FDTD method 
is proportional to the diameter of the computational cell, 
or A^i/'^ in d dimensions. With an ideal multigrid solver, 
FDFD can in principle be solved by 0(1) solver steps, but 
a simpler solver like conjugate gradient requires a number 
of steps proportional to the diameter as well [3(J . In both 
cases, the number of points to be solver on the surface S 
is 0{N^^^^'^). Hence, the overall complexity of the sim- 
plest implementations (not multigrid) is 0{N'^). We be- 
lieve that future boundary-element methods |30} 137] will 
achieve better efficiency, but such methods require con- 
siderable effort to implement and their implementation 
is specific to the homogeneous-medium Green's function, 
which depends on the boundary conditions, dimension- 
ality and types of materials considered f58j . 

Part II of this manuscript |52j . in preparation, will il- 
lustrate the method in various non-trivial two- and three- 
dimensional geometries, including dispersive dielectrics. 
In addition, we introduce an optimization of our method 
(based on a rapidly converging series expansion of the 
fields) that greatly speeds up the spatial integral of 
the stress tensor. We also compute forces in three- 
dimensional geometries with cylindrical symmetry, which 
allows us to take advantage of the cylindrical coordinates 
support in existing FDTD software [5T] and employ a 
two-dimensional computational cell. 



FDTD algorithms approximate both time and space 
by a discrete uniform mesh. Bearing aside the standard 
analysis of stability and convergence [48], this discretiza- 
tion will slightly modify the analysis in the preceding sec- 
tions. In particular, the use of a finite temporal grid (res- 
olution At) implies that all continuous time derivatives 
are now replaced by a finite-difference relation, which is 
most commonly taken to be a center difference: 



5/ Mr,t + At/2)~Mr,t-At/2) ^ (a) 
dt~ At - * ^ 



(35) 



where J{t) is an arbitrary function of time. The effect 
of temporal discretization is therefore to replace the lin- 
ear operator d/dt with d['^\ The representation of this 
operator is simple to compute in the frequency domain. 
Letting d'l'^'^ act on a Fourier component of of f{t) yields: 



where 



2 . f £.At\ „,5At 
— sm ^— e 2 
At V 2 / 



(36) 



(37) 



The effect of discretization on the system is thus to re- 
place by iS^d in the derivatives, which correspond to nu- 
merical dispersion arising from the ultraviolet (Nyquist) 
frequency cutoff it/ At. Note that ^ is still the frequency 
parameter governing the time dependence of the Fourier 
components of /(<) and ^ f in the limit of infinite 
resolution {At 0) . 

Because FDTD is convergent [£,d = £.+0{At^)], most of 
the analysis can be done (as in this paper) in the At ^ 
limit. However, care must be taken in computing g(t) 
because the Fourier transform of g{£,), Eq. (27 1, does not 
exist as At ^ 0. We must compute it in the finite At 
regime. In particular, the finite resolution requires, via 
Eq. (37), that we replace g{uj) in Eq. (27 1 by: 



5d(0 = 



(38) 



Note that the Jacobian factor dus/d^ involves cu and ^, 
not LUd and ^d, although of course the latter converges 
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to the former for At 0. The basic principle is that 
one must be careful to use the discrete analogues to con- 
tinuous solutions in cases where there is a divergence or 
regularization needed. This is the case for g{£,), but not 
for the Jacobian. 

Similarly, if one wished to subtract the vacuum Green's 
function from the Green's function, one needs to subtract 
the vacuum Green's function as computed in the dis- 
cretized vacuum. Such a subtraction is unnecessary if the 
stress tensor is integrated over a closed surface (vacuum 
contributions are constants that integrate to zero), but is 
useful in cases hke the parallel plates considered here. By 
subtracting the (discretized) vacuum Green's function, 
one can evaluate the stress tensor only for a single point 
between the plates, rather than for a "closed surface" 
with another point on the other side of the plates [5]. 



As was noted before Eq. (30 1, the Nyquist frequency 



n/At regularizes the frequency integrations, similar to 
other ultraviolet regularization schemes employed in 
Casimir force calculations ESI [60] . Because the total fre- 
quency integrand in Eq. M goes to zero for large ^ (due 
to cancellations occurring in the spatial integration and 
also due to the dissipation introduced in our approach), 
the precise nature of this regularization is irrelevant as 
long as At is sufficiently small (i.e., at high enough reso- 
lution) . 



B. The Magnetic Correlation Function 

One way to compute the magnetic correlation func- 
tion is by taking spatial derivatives of the electric Green's 
function by Eq. (|8]), but this is inefficient because a nu- 
merical derivative involves evaluating the electric Green's 
function at multiple points. Instead, we compute the 
magnetic Green's function directly, finding the magnetic 
field in response to a magnetic current. This formulation, 
however, necessitates a change in the choice of vector po- 
tentials Eqs. (|3||4| as well as a switch from an electric 
to magnetic conductivity, for reasons explained in this 
section. 



Equations (|3||4|) express the magnetic field B as the 
curl of the vector potential A^, enforcing the constraint 
that B is divergence-free (no magnetic charge). However, 
this is no longer true when there is a magnetic current, 
as can be seen by taking the divergence of both sides 
of Faraday's law with a magnetic current J, d'B/dt = 
— V X E — J, since V • J 7^ for a point-dipole current J. 
Instead, since there need not be any free electric charge in 
the absence of an electric current source, one can switch 
to a new vector potential A^ such that 



e£;,(r,c.) = (Vx)^^.^f(r,^) 



(39) 
(40) 



The desired correlation function is then given, analogous 
to Eq. 0, by 

{K,{r,u;)H,{r',u;)) - ^c.^ Im G^(c^; r, r'), (41) 

where the photon magnetic Green's function solves 
[similar to Eq. ^] 



V X -y-^ — ^ V X — cj^u(r, uj) 
e(r,w) 



Gf(c.;r,r') = <5(r-r')e,. 

(42) 



Now, all that remains is to map Eq. (42 1 onto an equiv- 



alent real-frequency (^) system that can be evaluated in 
the time domain, similar to Sec. II C for uj(£^) given by 
Eq. (18). There are at least two ways to accomplish 
this. One possibility, which we have adopted in this 
paper, is to define an effective magnetic permeability 
/ic = /^"-^^(O/C^i corresponding to a magnetic conduc- 
tivity, similar to Eq. (13 1. Combined with Eq. (18 1, this 



directly yields a magnetic conductivity as in Eq. ( 20 ) . 

A second possibility is to divide both sides of Eq. ( 42 1 
by w^/^^ = 1 + zcr/^, and absorb the 1 + ia/S, factor into 
£ via Eq. ( [13^ . That is, one can compute the magnetic 
correlation function via the magnetic field in response to 
a magnetic current with an electric conductivity. How- 
ever, the magnetic current in this case has a frequency 
response that is divided by 1 -I- icr/C, which is simply 
a rescaling of J^{£,) in Eq. (22). There is no particular 



computational advantage to this alternative, but for an 
experimental realization |50j . an electric conductivity is 
considerably more attractive. [Note that rescaling J'(f) 
by l + ia/uj will yield a new g{^) in Eq. (27), correspond- 
ing to a new g(t) that exhibits slower decay.] 



C. Material Dispersion 

In this section, we extend the time-domain formalism 
presented above to cases where the dielectric permittivity 
of the medium of interest is dispersive. To begin with, 
note that in this case the dissipative, complex dielectric 
Ec of Eq. (43 1 is given by: 



£c(r,0 



(43) 



where e(r, uj{£,)) denotes the permittivity of the geometry 
of interest evaluated over the complex contour lu{S,)- 

This complex dielectric manifests itself as a convolu- 
tion in the time-domain equations of motion, i.e. in gen- 
eral, D(t) = J dt'edt - t')E{t'). The standard way to 
implement this in FDTD is to employ an auxiliary equa- 
tion of motion for the polarization [57 . For the particular 
contour chosen in this paper [Eq. ([T8|], the conductivity 
term already includes the prefactor uP' and therefore 
one need only add the dispersion due to e(r,w(f)). 

The only other modification to the method comes from 
the dependence of r^(^) in Eq. (25) on e. We remind 
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the reader that our definition of T was motivated by our 



desire to interpret Eq. (24 1 as the Fourier transform of 



the convolution of two quantities, and thus to express the 
Casimir force directly in terms of the electric and magen- 
tic fields E(t) and H(i), respectively. A straightforward 
generalization of Eq. p5| ) to dispersive media entails set- 
ting £(r) e{r,uj). However, in this case, the Fourier 



transform of Eq. ( 25 ) would be given by a convolution 
of E(^) and e(r g S,uj{£^)) in the time domain, making 
it impossible to obtain T^{t) directly in terms of E(i). 
This is not a problem however, because the stress ten- 
sor must be evaluated over a surface S that lies entirely 



within a uniform medium (otherwise, S would cross a 
boundary and interpreting the result as a force on par- 
ticular objects inside S would be problematic). The di- 
electric appearing in Eq. ( 25 1 is then at most a function 
of cj(^), i.e. e(r g S,uj) — e{uj), which imphcs that we 
can simply absorb this factor into g{S,), modifying the 
numerical integral of Eq. (30 1. Furthermore, the most 



common case considered in Casimir-force calculations is 
one in which the stress tensor is evaluated in vacuum, 
i.e. e(r G 5, w) = 1, and thus dispersion does not modify 
g{^) at all. 
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